<?php
//A library of Stats functions.  Version 1.8, June 20, 2007

global $allowedmacros;
array_push($allowedmacros,"nCr","nPr","mean","stdev","percentile","quartile","median","freqdist","frequency","histogram","fdhistogram","fdbargraph","normrand","boxplot","normalcdf","tcdf","invnormalcdf","invtcdf","invtcdf2","linreg","countif","binomialpdf","binomialcdf","chi2cdf","invchi2cdf","fcdf","invfcdf");

//nCr(n,r)
//The Choose function
function nCr($n,$r){
   if ($r > $n)
     return false;
   if (($n-$r) < $r)
     return nCr($n,($n-$r));
   $return = 1;
   for ($i=0;$i < $r;$i++){
     $return *= ($n-$i)/($i+1);
   }
   return $return;
}

//nPr(n,r)
//The Permutations function
function nPr($n,$r){
   if ($r > $n)
     return false;
   if ($r==0) 
     return 1;
   $return = 1;
   $i = $n-$r;
   while ($i<$n) {
     $return *= ++$i;
   }
   return $return;
}

//mean(array)
//Finds the mean of an array of numbers
function mean($a) {
	return (array_sum($a)/count($a));	
}

//variance(array)
//the (sample) variance of an array of numbers
function variance($a) {
	$v = 0;
	$mean = mean($a);
	foreach ($a as $x) {
		$v += pow($x-$mean,2);
	}
	return ($v/(count($a)-1));
}

//stdev(array)
//the (sample) standard deviation of an array of numbers
function stdev($a) {
	return sqrt(variance($a));
}

//percentile(array,percentile)
//example: percentile($a,30) would find the 30th percentile of the data
//method based on Triola
function percentile($a,$p) {
	sort($a, SORT_NUMERIC);
	if ($p==0) {
		return $a[0];
	} else if ($p==100) {
		return $a[count($a)-1];
	}

	$l = $p*count($a)/100;
	if (floor($l)==$l) {
		return (($a[$l-1]+$a[$l])/2);
	} else {
		return ($a[ceil($l)-1]);
	}
}

//quartile(array,quartile)
//finds the 0 (min), 1st, 2nd (median), 3rd, or 4th (max) quartile of an
//array of numbers.  Calculates using percentiles.
function quartile($a,$q) {
	return percentile($a,$q*25);	
}

//median(array)
//returns the median of an array of numbers
function median($a) {
	return percentile($a,50);	
}

//freqdist(array,label,start,classwidth)
//display macro.  Returns an HTML table that is a frequency distribution of
//the data
// array: array of data values
// label: name of data values
// start: first lower class limit
// classwidth: width of the classes
function freqdist($a,$label,$start,$cw) {
	sort($a, SORT_NUMERIC);
	$x = $start;
	$curr = 0;
	$out = "<table class=stats><thead><tr><th>$label</th><th>Freq</th></tr></thead>\n<tbody>\n";
	while ($x <= $a[count($a)-1]) {
		$out .= "<tr><td>`$x <= x < ";
		$x += $cw;
		$out .= "$x`</td><td>";
		$i = $curr;
		while (($a[$i] < $x) && ($i < count($a))) {
			$i++;
		}
		$out .= ($i-$curr) . "</td></tr>\n";
		$curr = $i;
	}
	$out .= "</tbody></table>\n";
	return $out;
}

//frequency(array,start,classwidth)
//Returns an array of frequencies for the data grouped into classes
// array: array of data values
// start: first lower class limit
// classwidth: width of the classes
function frequency($a,$start,$cw) {
	sort($a, SORT_NUMERIC);
	$x = $start;
	$curr = 0;
	while ($x <= $a[count($a)-1]) {
		$x += $cw;
		$i = $curr;
		while (($a[$i] < $x) && ($i < count($a))) {
			$i++;
		}
		$out[] = ($i-$curr);
		$curr = $i;
	}
	return $out;
}

//countif(array,condition)
//Returns count of items in array that meet condition
// array: array of data values
// condition: a condition, using x for data values
//Example: countif($a,"x<3 && x>2")
function countif($a,$ifcond) {
	$rsnoquote = preg_replace('/"[^"]*"/','""',$ifcond);
	$rsnoquote = preg_replace('/\'[^\']*\'/','\'\'',$rsnoquote);
	if (preg_match_all('/([$\w]+)\s*\([^\)]*\)/',$rsnoquote,$funcs)) {
		$ismath = true;
		for ($i=0;$i<count($funcs[1]);$i++) {
			if (strpos($funcs[1][$i],"$")===false) {
				if (!in_array($funcs[1][$i],$allowedmacros)) {
					echo "{$funcs[1][$i]} is not an allowed function<BR>\n";
					return false;
				}
			}
		}
	}
	$ifcond = str_replace('!=','#=',$ifcond);
	$ifcond = mathphp($ifcond,'x');
	$ifcond = str_replace('#=','!=',$ifcond);
	$ifcond = str_replace('(x)','($x)',$ifcond);
	$iffunc = create_function('$x','return('.$ifcond.');');
	
	$cnt = 0;
	foreach ($a as $v) {
		if ($iffunc($v)) {
			$cnt++;
		}
	}
	return $cnt;
}

//histogram(array,label,start,classwidth,[labelstart,upper,width,height])
//display macro.  Creates a histogram from a data set
// array: array of data values
// label: name of data values
// start: first lower class limit
// classwidth: width of the classes
// labelstart (optional): value to start axis labeling at.  Defaults to start
// upper (optional): first upper class limit.  Defaults to start+classwidth
// width,height (optional): width and height in pixels of graph
function histogram($a,$label,$start,$cw,$startlabel=false,$upper=false,$width=300,$height=200) {
	sort($a, SORT_NUMERIC);
	$x = $start;
	$curr = 0;
	$alt = "Histogram for $label <table class=stats><thead><tr><th>Label on left of box</th><th>Frequency</th></tr></thead>\n<tbody>\n";
	$maxfreq = 0;
	if ($upper===false) {
		$dx = $cw;
		$dxdiff = 0;
	} else {
		$dx = $upper - $start;
		$dxdiff = $cw-$dx;
	}

	while ($x <= $a[count($a)-1]) {
		$alt .= "<tr><td>$x</td>";
		$st .= "rect([$x,0],";
		$x += $dx;
		$st .= "[$x,";
		$i = $curr;
		while (($a[$i] < $x) && ($i < count($a))) {
			$i++;
		}
		if (($i-$curr)>$maxfreq) { $maxfreq = $i-$curr;}
		$alt .= "<td>" . ($i-$curr) . "</td></tr>\n";
		$st .= ($i-$curr) . "]);";
		$curr = $i;
		$x += $dxdiff;
	}
	$alt .= "</tbody></table>\n";
	if ($GLOBALS['sessiondata']['graphdisp']==0) {
		return $alt;
	}
	$outst = "setBorder(10);  initPicture(". ($start-.1*($x-$start)) .",$x,". (-.1*$maxfreq) .",$maxfreq);";
	if ($maxfreq>100) {$step = 20;} else if ($maxfreq > 50) { $step = 10; } else if ($maxfreq > 20) { $step = 5;} else if ($maxfreq>10) { $step = 2; } else {$step=1;}
	if ($startlabel===false) {
		$outst .= "axes($cw,$step,1,1000,$step); fill=\"blue\"; text([". ($start + .5*($x-$start))  .",". (-.1*$maxfreq) . "],\"$label\");";
	} else {
		$outst .= "axes(1000,$step,1,1000,$step); fill=\"blue\"; text([". ($start + .5*($x-$start))  .",". (-.1*$maxfreq) . "],\"$label\");";
		$x = $startlabel;
		$tm = -.02*$maxfreq;
		$tx = .02*$maxfreq;
		while ($x <= $a[count($a)-1]) {
			$outst .= "line([$x,$tm],[$x,$tx]); text([$x,0],\"$x\",\"below\");";
			$x+= $cw;
		}

	}
	$outst .= $st;
	return showasciisvg($outst,$width,$height);
}

//fdhistogram(freqarray,label,start,cw,[labelstart,upper,width,height])
//display macro.  Creates a histogram from frequency array
// freqarray: array of frequencies
// label: name of data values
// start: first lower class limit
// classwidth: width of the classes
// labelstart (optional): value to start axis labeling at.  Defaults to start
// upper (optional): first upper class limit.  Defaults to start+classwidth
// width,height (optional): width and height in pixels of graph
function fdhistogram($freq,$label,$start,$cw,$startlabel=false,$upper=false,$width=300,$height=200) {
	$x = $start;
	$alt = "Histogram for $label <table class=stats><thead><tr><th>Label on left of box</th><th>Frequency</th></tr></thead>\n<tbody>\n";
	$maxfreq = 0;
	if ($upper===false) {
		$dx = $cw;
		$dxdiff = 0;
	} else {
		$dx = $upper - $start;
		$dxdiff = $cw-$dx;
	}
	for ($curr=0; $curr<count($freq); $curr++) {
		$alt .= "<tr><td>$x</td><td>{$freq[$curr]}</td></tr>";
		$st .= "rect([$x,0],";
		$x += $dx;
		$st .= "[$x,{$freq[$curr]}]);";
		if ($freq[$curr]>$maxfreq) { $maxfreq = $freq[$curr];}
		$x += $dxdiff;
	}
	$alt .= "</tbody></table>\n";
	if ($GLOBALS['sessiondata']['graphdisp']==0) {
		return $alt;
	}
	$outst = "setBorder(10);  initPicture(". ($start-.1*($x-$start)) .",$x,". (-.1*$maxfreq) .",$maxfreq);";
	if ($maxfreq>100) {$step = 20;} else if ($maxfreq > 50) { $step = 10; } else if ($maxfreq > 20) { $step = 5;} else {$step=1;}
	if ($startlabel===false) {
		$outst .= "axes($cw,$step,1,1000,$step); fill=\"blue\"; text([". ($start + .5*($x-$start))  .",". (-.1*$maxfreq) . "],\"$label\");";
	} else {
		$outst .= "axes(1000,$step,1,1000,$step); fill=\"blue\"; text([". ($start + .5*($x-$start))  .",". (-.1*$maxfreq) . "],\"$label\");";
		$x = $startlabel;
		$tm = -.02*$maxfreq;
		$tx = .02*$maxfreq;
		for ($curr=0; $curr<count($freq); $curr++) {
			$outst .= "line([$x,$tm],[$x,$tx]); text([$x,0],\"$x\",\"below\");";
			$x+=$cw;
		}
	}
	//$outst .= "axes($cw,$step,1,1000,$step); fill=\"blue\"; text([". ($start + .5*($x-$start))  .",". (-.1*$maxfreq) . "],\"$label\");";
	$outst .= $st;
	return showasciisvg($outst,$width,$height);
}

//fdbargraph(barlabels,freqarray,label,[width,height])
//barlabels: array of labels for the bars
//freqarray: array of frequencies/heights for the bars
//label: general label for bars
//width,height (optional): width and height for graph
function fdbargraph($bl,$freq,$label,$width=300,$height=200) {
	$alt = "Histogram for $label <table class=stats><thead><tr><th>Bar Label</th><th>Frequency/Height</th></tr></thead>\n<tbody>\n";
	$start = 0;
	$x = $start+1;
	$maxfreq = 0;
	for ($curr=0; $curr<count($bl); $curr++) {
		$alt .= "<tr><td>{$bl[$curr]}</td><td>{$freq[$curr]}</td></tr>";
		$st .= "rect([$x,0],";
		$x += 2;
		$st .= "[$x,{$freq[$curr]}]);";
		$x -= 1;
		$st .= "text([$x,0],\"{$bl[$curr]}\",\"below\");";
		$x += 1;
		if ($freq[$curr]>$maxfreq) { $maxfreq = $freq[$curr];}
	}
	$alt .= "</tbody></table>\n";
	if ($GLOBALS['sessiondata']['graphdisp']==0) {
		return $alt;
	}
	$x++;
	$outst = "setBorder(10);  initPicture(". ($start-.1*($x-$start)) .",$x,". (-.1*$maxfreq) .",$maxfreq);";
	if ($maxfreq>100) {$step = 20;} else if ($maxfreq > 50) { $step = 10; } else if ($maxfreq > 20) { $step = 5;} else {$step=1;}
	$outst .= "axes(1000,$step,1,1000,$step); fill=\"blue\"; text([". ($start + .5*($x-$start))  .",". (-.12*$maxfreq) . "],\"$label\");";
	
	//$outst .= "axes($cw,$step,1,1000,$step); fill=\"blue\"; text([". ($start + .5*($x-$start))  .",". (-.1*$maxfreq) . "],\"$label\");";
	$outst .= $st;
	return showasciisvg($outst,$width,$height);
}

//normrand(mu,sigma,n)
//returns an array of n random numbers that are normally distributed with given
//mean mu and standard deviation sigma.  Uses the Box-Muller transform.
function normrand($mu,$sig,$n) {
	for ($i=0;$i<ceil($n/2);$i++) {
		do {
			$a = rand(-32768,32768)/32768;
			$b = rand(-32768,32768)/32768;
			$r = $a*$a+$b*$b;
			$count++;
		} while ($r==0||$r>1);
		$r = sqrt(-2*log($r)/$r); 
		$z[] = $sig*$a*$r + $mu;
		$z[] = $sig*$b*$r + $mu;
	}
	if ($n%2==0) {
		return $z;
	} else {
		return (array_slice($z,0,count($z)-1));
	}
}

//boxplot(array,axislabel,[datalabel])
//draws a boxplot based on the data in array, with given axislabel
//and optionally a datalabel (to topleft of boxplot)
//array and datalabel can also be an array of dataarrays and
//array of datalabels to do comparative boxplots
function boxplot($arr,$label) {
	if (func_num_args()>2) {
		$dlbls = func_get_arg(2);
	}
	if (is_array($arr[0])) { $multi = count($arr);} else {$multi = 1;}
	$st = '';
	$bigmax = -10000000;
	$bigmin = 100000000;
	$alt = "Boxplot, axis label: $label. ";
	for ($i=0;$i<$multi;$i++) {
		if ($multi==1) { $a = $arr;} else {$a = $arr[$i];}
		sort($a,SORT_NUMERIC);
		$max = $a[count($a)-1];
		if ($max>$bigmax) { $bigmax = $max;}
		if ($a[0]<$bigmin) {$bigmin = $a[0];}
		$q1 = percentile($a,25);
		$q2 = percentile($a,50);
		$q3 = percentile($a,75);
		$yl = 2+5*$i;
		$ym = $yl+1;
		$yh = $ym+1;
		$st .= "line([$a[0],$yl],[$a[0],$yh]); rect([$q1,$yl],[$q3,$yh]); line([$q2,$yl],[$q2,$yh]);";
		$st .= "line([$max,$yl],[$max,$yh]); line([$a[0],$ym],[$q1,$ym]); line([$q3,$ym],[$max,$ym]);";
		$alt .= 'Boxplot ';
		if (isset($dlbls[$i])) {
			$alt .= "for {$dlbls[$i]}";
		} else {
			$alt .= ($i+1);
		}
		$alt .= ": Left whisker at {$a[0]}. Leftside of box at $q1. Line through box at $q2.  Rightside of box at $q3. Right whisker at $max.\n"; 
		
	}
	if ($GLOBALS['sessiondata']['graphdisp']==0) {
		return $alt;
	}
	$outst = "setBorder(15); initPicture($bigmin,$bigmax,-3,".(5*$multi).");";
	$dw = $bigmax-$bigmin;

	if ($dw>100) {$step = 20;} else if ($dw > 50) { $step = 10; } else if ($dw > 20) { $step = 5;} else {$step=1;}
	$outst .= "axes($step,100,1);";
	$outst .= "text([". ($bigmin+.5*$dw) . ",-3],\"$label\");";
	if (isset($dlbls)) {
		for ($i=0;$i<$multi;$i++) {
			if ($multi>1) { $dlbl = $dlbls[$i];} else {$dlbl = $dlbls;}
			$st .= "text([$bigmin,". (5+5*$i)  ."],\"$dlbl\",\"right\");";
		}
	}
	$outst .= $st;
	return showasciisvg($outst,400,50+50*$multi);
}

//normalcdf(z,[dec])
//calculates the area under the standard normal distribution to the left of the
//z-value z, to dec decimals (defaults to 4) 
//based on someone else's code - can't remember whose!
function normalcdf($ztest,$dec=4) {
	$eps = pow(.1,$dec);
	$eps2 = pow(.1,$dec+2);
	
	$ds = 1;
	$s = 0;
	$i = 0;
	$z = abs($ztest);
	$fact = 1;
	while (abs($ds)>$eps2) {
		$ds = pow(-1,$i)*pow($z,2.0*$i+1.0)/(pow(2.0,$i)*$fact*(2.0*$i+1.0));
		$s += $ds;
		$i++;
		$fact *= $i;
		if (abs($s)<$eps) {
			break;
		}
	}
	
	$s *= 0.39894228;
	$s = round($s,$dec);
	if ($ztest > 0) {
		$pval = .5 + $s;
	} else {
		$pval = .5 - $s;
	}
	if ($pval < $eps) {
		$pval = $eps;
	} else if ($pval > (1-$eps)) {
		$pval = round(1-$eps,$dec);
	}
	return $pval;
}

//tcdf(t,df,[dec])
//calculates the area under the t-distribution with "df" degrees of freedom
//to the left of the t-value t
//based on someone else's code - can't remember whose!
function tcdf($ttest,$df,$dec=4) {
	$eps = pow(.1,$dec);
	$eps2 = pow(.1,$dec+2);
	
	$t = abs($ttest);
	if ($df > 0) {
	$k3 = 0;
	$c3 = 0;
	$a3 = $t/sqrt($df);
	$b3 = 1+pow($a3,2);
	$y = 0.5;
	if (abs(floor($df/2)*2 - $df) < $eps) { 
		$c3 = $a3/(2*pow($b3,0.5));
		$k3 = 0;
	} else {
		$c3 = $a3/($b3*M_PI);
		$y += atan($a3)/M_PI;
		$k3 = 1;
	}
	if ($df > 1) {
		$k3 += 2;
		$y += $c3;
		$c3 *= ($k3-1)/($k3*$b3);
		while ($k3 < $df) {
			$k3 += 2;
			$y += $c3;
			$c3 *= ($k3-1)/($k3*$b3);
		}
	}
	if ($ttest > 0) {
		$pval = round($y,$dec);
	} else {
		$pval = round(1-$y,$dec);
	}


	if ($pval > (1-$eps)) {
		$pval = round(1-$eps,$dec);
	} else if ($pval < $eps) {
		$pval = $eps;
	}
	return $pval;
	} else {return false;}
}

//invnormalcdf(p,[dec])
//Inverse Normal CDF
//finds the z-value with a left-tail area of p, to dec decimals (default 5)
// from Odeh & Evans. 1974. AS 70. Applied Statistics. 23: 96-97
function invnormalcdf($p,$dec=5) {
   
      $p0 = -0.322232431088;
      $p1 = -1.0;
      $p2 = -0.342242088547;
      $p3 = -0.0204231210245;
      $p4 = -0.453642210148e-4;
      $q0 =  0.0993484626060;
      $q1 =  0.588581570495;
      $q2 =  0.531103462366;
      $q3 =  0.103537752850;
      $q4 =  0.38560700634e-2;
      
   if ($p < 0.5) { $pp = $p; }  else  {$pp = 1 - $p; }


   if ($pp < 1E-12) {
      $xp = 99;
   } else {
      $y = sqrt(log(1/($pp*$pp)));
      $xp = $y + (((($y * $p4 + $p3) * $y + $p2) * $y + $p1) * $y + $p0) / (((($y * $q4 + $q3) * $y + $q2) * $y + $q1) * $y + $q0);
   }
   $xp = round($xp,$dec);
   if ($p < 0.5) { return (-1*$xp); }  else { return  $xp; }
}
   
//invtcdf(p,df,[dec])
//the inverse Student's t-distribution 
//computes the t-value with a left-tail probability of p, with df degrees of freedom
//to dec decimal places (default 4)
// from Algorithm 396: Student's t-quantiles by G.W. Hill  Comm. A.C.M., vol.13(10), 619-620, October 1970
function invtcdf($p,$ndf,$dec=4) {
	$half_pi = M_PI/2;
	$eps = 1e-12;
	$origp = $p;
	if ($p > 0.5)  {$p = 1 - $p; }
	$p = $p*2;
	
	if ($ndf < 1 || $p > 1 || $p <= 0) {
	  echo "error in params";
	  return false;
	}
	
	if ( abs($ndf - 2) < $eps ) {
	  $fn_val = SQRT(2 / ( $p * (2 - $p) ) - 2);
	  if ($origp<.5) {return (-1*$fn_val);} else { return $fn_val;}
	} else if ($ndf < 1+$eps) {
	  $prob = $p * $half_pi;
	  $fn_val = cos($prob) / sin($prob);
	  if ($origp<.5) {return (-1*$fn_val);} else { return $fn_val;}
	} else {
	  $a = 1/ ($ndf - 0.5);
	  $b = 48/ pow($a,2);
	  $c = ((20700 * $a / $b - 98) * $a - 16) * $a + 96.36;
	  $d = ((94.5 / ($b + $c) - 3) / $b + 1) * sqrt($a * $half_pi)* $ndf;
	  $x = $d * $p;
	  $y = pow($x,(2/ $ndf));
	  
	  if ($y > 0.05 + $a) {
	    $x = invnormalcdf($origp);
	    $y = pow($x,2);
	    if ($ndf < 5) { $c = $c + 0.3 * ($ndf - 4.5) * ($x + 0.6);}
	    $c = (((0.05 * $d * $x - 5) * $x - 7) * $x - 2) * $x + $b+ $c;

	    $y = (((((0.4*$y + 6.3) * $y + 36) * $y + 94.5) / $c - $y - 3) / $b + 1) * $x;
	    $y = $a * pow($y,2);
	    if ($y > 0.002) {
	      $y = exp($y) - 1;
	    } else {
	      $y = 0.5 * pow($y,2) + $y;
	    }
	  } else {
	    
	    $y = ((1 / ((($ndf + 6) / ($ndf * $y) - 0.089 * $d - 0.822) * ($ndf + 2) * 3) + 0.5 / ($ndf + 4)) * $y - 1) * ($ndf + 1) / ($ndf + 2) + 1 / $y;
	  }
	}
	$fn_val = round( sqrt($ndf * $y) , $dec);
	if ($origp<.5) {return (-1*$fn_val);} else { return $fn_val;}
}

//linreg(xarray,yarray)
//Computes the linear correlation coefficient, and slope and intercept of
//regression line, based on array/list of x-values and array/list of y-values
//Returns as array:  r,slope,intercept
function linreg($xarr,$yarr) {
	if (!is_array($xarr)) { $xarr = explode(',',$xarr);}	
	if (!is_array($yarr)) { $yarr = explode(',',$yarr);}
	if (count($xarr)!=count($yarr)) { 
		echo "Error: linreg requires xarray length = yarray length";
		return false;
	}
	$sx = array_sum($xarr);
	$sy = array_sum($yarr);
	$sxx=0; $syy=0; $sxy = 0;
	for ($i=0; $i<count($xarr); $i++) {
		$sxx += $xarr[$i]*$xarr[$i];
		$syy += $yarr[$i]*$yarr[$i];
		$sxy += $xarr[$i]*$yarr[$i];
	}
	$n = count($xarr);
	$r = ($n*$sxy - $sx*$sy)/(sqrt($n*$sxx-$sx*$sx)*sqrt($n*$syy-$sy*$sy));
	$m = ($n*$sxy - $sx*$sy)/($n*$sxx - $sx*$sx);
	$b = ($sy - $sx*$m)/$n;
	return array($r,$m,$b);
}

//binomialpdf(N,p,x)
//Computes the probability of x successes out of N trials
//where each trial has probability p of success
function binomialpdf($N,$p,$x) {
	return (nCr($N,$x)*pow($p,$x)*pow(1-$p,$N-$x));
}


//binomialcdf(N,p,x)
//Computes the probably of &lt;=x successes out of N trials
//where each trial has probability p of success
function binomialcdf($N,$p,$x) {
	$out = 0;
	for ($i=0;$i<=$x;$i++) {
		$out += binomialpdf($N,$p,$i);
	}
	return $out;
}

//chi2cdf(x,df)
//Computes the area to the left of x under the chi-squared disribution
//with df degrees of freedom
function chi2cdf($x,$a) {
	return gamma_cdf(0.5*$x,0.0,1.0,0.5*$a);
}

//invchi2cdf(p,df)
//Compuates the x value with left-tail probability p under the 
//chi-squared distribution with df degrees of freedom
function invchi2cdf($cdf,$a) {
	$aa = 0.6931471806;
  $c1 = 0.01;
  $c2 = 0.222222;
  $c3 = 0.32;
  $c4 = 0.4;
  $c5 = 1.24;
  $c6 = 2.2;
  $c7 = 4.67;
  $c8 = 6.66;
  $c9 = 6.73;
  $c10 = 13.32;
  $c11 = 60.0;
  $c12 = 70.0;
  $c13 = 84.0;
  $c14 = 105.0;
  $c15 = 120.0;
  $c16 = 127.0;
  $c17 = 140.0;
  $c18 = 175.0;
  $c19 = 210.0;
  $c20 = 252.0;
  $c21 = 264.0;
  $c22 = 294.0;
  $c23 = 346.0;
  $c24 = 420.0;
  $c25 = 462.0;
  $c26 = 606.0;
  $c27 = 672.0;
  $c28 = 707.0;
  $c29 = 735.0;
  $c30 = 889.0;
  $c31 = 932.0;
  $c32 = 966.0;
  $c33 = 1141.0;
  $c34 = 1182.0;
  $c35 = 1278.0;
  $c36 = 1740.0;
  $c37 = 2520.0;
  $c38 = 5040.0;
  $cdf_max = 0.999998;
  $cdf_min = 0.000002;
  $e = 0.0000005;

  $it_max = 20;
  if ($cdf < $cdf_min) {
	  echo "p < min - can't do it!";
	  return 0;
  }
  if ($cdf_max < $cdf) {
	  echo "p > max - can't do it";
	  return 0;
  }
  $xx = 0.5*$a;
  $c = $xx-1.0;
  $g = gamma_log($a/2.0);
  if ($a < -$c5*log($cdf)) {
	  $ch = pow(($cdf*$xx*exp($g+$xx*$aa)),(1.0/$xx));
	  if ($ch < $e) {
		  $x = $ch;
		  return $x;
	  }
  } else if ($a <= $c3) {
	  $ch = $c4;
	  $a2 = log(1.0 - $cdf);
	  while (1) {
		  $q = $ch;
		  $p1 = 1.0+$ch*($c7 + $ch);
		  $p2 = $ch*($c9 + $ch*($c8+$ch));
		  $t = -0.5+($c7 + 2.0*$ch)/$p1 - ($c9 + $ch*($c10 + 3.0*$ch)) / $p2;
		  $ch = $ch - (1.0 - exp($a2 + $g + 0.5*$ch + $c*$aa)*$p2/$p1)/$t;
		  if (abs($q/$ch - 1.0) <= $c1) {
			  break;
		  }
	  }
  } else {
	  $x2 = invnormalcdf($cdf);
	  $p1 = $c2/$a;
	  $ch = $a*pow(($x2*sqrt($p1) + 1.0 - $p1),3);
	  if ($c6*$a + 6.0 < $ch) {
		  $ch = -2.0*(log(1.0 - $cdf) - $c*log(0.5*$ch)+$g);
	  }
  }
  for ($i=1;$i<=$it_max;$i++) {
	  $q = $ch;
	  $p1 = 0.5*$ch;
	  $p2 = $cdf - gamma_inc($xx,$p1);
	  $t = $p2*exp($xx*$aa + $g + $p1 - $c*log($ch));
	  $b = $t/$ch;
	  $a2 = 0.5*$t - $b*$c;
	  $s1 = ($c19 + $a2 * ($c17 + $a2 * ($c14 + $a2 *($c13 + $a2 * ($c12 + $a2 * $c11))))) / $c24;
	  $s2 = ($c24 + $a2 * ($c29 + $a2 * ($c32 + $a2 *($c33 + $a2 * $c35)))) / $c37;
	  $s3 = ($c19 + $a2 * ($c25 + $a2 * ($c28 + $a2 * $c31)))/$c37;
	  $s4 = ($c20 + $a2 * ($c27 + $a2 * $c34) + $c * ($c22 + $a2 * ($c30 + $a2 *$c36)))/$c38;
	  $s5 = ($c13 + $c21 + $a2 + $c * ($c18 + $c26 * $a2))/$c37;
	  $s6 = ($c16 + $c*($c23 + $c16*$c))/$c38;
	  $ch = $ch + $t*(1.0 + 0.5*$t*$s1 - $b*$c * ($s1-$b*($s2-$b*($s3-$b*($s4-$b*($s5-$b*$s6))))));
	  if ($e < abs($q/$ch - 1.0)) {
		  $x = $ch;
		  return ($x);
	  }
  }
  $x = $ch;
  echo "Convergence not reached in invchisquare";
  return $x;
}

function gamma_cdf($x,$a,$b,$c) {
	return gamma_inc($c,($x-$a)/$b);
}
function gamma_inc($p,$x,$dec=4) {
	$exp_arg_min = -88.0;
	$overflow = 1e37;
	$plimit = 1000;
	$tol = 1e-7;
	$xbig = 1e8;
	$value = 0.0;
	
	if ($plimit<$p) {
		$pn1 = 3.0*sqrt($p)*(pow($x/$p,1.0/3.0) + 1.0/(9.0*$p)-1.0);
		$value = normalcdf($pn1,$dec);
		return $value;
	}
	if ($xbig<$x) {
		return 1.0;
	}
	if ($x<=1.0 || $x<$p) {
		$arg = $p*log($x) - $x - gamma_log($p+1.0);
		$c = 1.0;
		$value = 1.0;
		$a = $p;
		while ($c>$tol) {
			$a += 1.0;
			$c = $c*$x/$a;
			$value += $c;
		}
		$arg += log($value);
		if ($exp_arg_min <= $arg) {
			$value = exp($arg);
		} else {
			$value = 0.0;
		}
	} else {
		$arg = $p*log($x)-$x-gamma_log($p);
		$a = 1.0 - $p;
		$b = $a+$x+1.0;
		$c = 0.0;
		$pn1 = 1.0;
		$pn2 = $x;
		$pn3 = $x+1.0;
		$pn4 = $x*$b;
		$value = $pn3/$pn4;
		while (1) {
			$a += 1.0;
			$b += 2.0;
			$c += 1.0;
			$pn5 = $b*$pn3 - $a*$c*$pn1;
			$pn6 = $b*$pn4 - $a*$c*$pn2;
			
			if (0 < abs($pn6)) {
				$rn = $pn5/$pn6;
				if (abs($value-$rn) <= min($tol,$tol*$rn)) {
					$arg += log($value);
					if ($exp_arg_min <= $arg) {
						$value = 1.0 - exp($arg);
					} else {
						$value = 1.0;
					}
					return $value;
				}
				$value = $rn;
			}
			$pn1 = $pn3;
			$pn2 = $pn4;
			$pn3 = $pn5;
			$pn4 = $pn6;
			if ($overflow <= abs($pn5)) {
				$pn1 = $pn1/$overflow;
				$pn2 = $pn2/$overflow;
				$pn3 = $pn3/$overflow;
				$pn4 = $pn4/$overflow;
			}
		}
	}
	return $value;
}

function gamma_log($x) {
 $c = array(
    -1.910444077728E-03, 
     8.4171387781295E-04, 
    -5.952379913043012E-04, 
     7.93650793500350248E-04, 
    -2.777777777777681622553E-03, 
     8.333333333333333331554247E-02, 
     5.7083835261E-03 );
  $d1 = - 5.772156649015328605195174E-01;
  $d2 =   4.227843350984671393993777E-01;
  $d4 =   1.791759469228055000094023;
  $frtbig = 1.42E+09;
  $p1 = array(
    4.945235359296727046734888, 
    2.018112620856775083915565E+02, 
    2.290838373831346393026739E+03, 
    1.131967205903380828685045E+04, 
    2.855724635671635335736389E+04, 
    3.848496228443793359990269E+04, 
    2.637748787624195437963534E+04, 
    7.225813979700288197698961E+03 );
  $p2 = array(
    4.974607845568932035012064, 
    5.424138599891070494101986E+02, 
    1.550693864978364947665077E+04, 
    1.847932904445632425417223E+05, 
    1.088204769468828767498470E+06, 
    3.338152967987029735917223E+06, 
    5.106661678927352456275255E+06, 
    3.074109054850539556250927E+06 );
   $p4 = array(
    1.474502166059939948905062E+04, 
    2.426813369486704502836312E+06, 
    1.214755574045093227939592E+08, 
    2.663432449630976949898078E+09, 
    2.940378956634553899906876E+010,
    1.702665737765398868392998E+011,
    4.926125793377430887588120E+011, 
    5.606251856223951465078242E+011 );
  $pnt68 = 0.6796875;
  $q1 = array(
    6.748212550303777196073036E+01, 
    1.113332393857199323513008E+03, 
    7.738757056935398733233834E+03, 
    2.763987074403340708898585E+04, 
    5.499310206226157329794414E+04, 
    6.161122180066002127833352E+04, 
    3.635127591501940507276287E+04, 
    8.785536302431013170870835E+03 );
  $q2 = array(
    1.830328399370592604055942E+02, 
    7.765049321445005871323047E+03, 
    1.331903827966074194402448E+05, 
    1.136705821321969608938755E+06, 
    5.267964117437946917577538E+06, 
    1.346701454311101692290052E+07, 
    1.782736530353274213975932E+07, 
    9.533095591844353613395747E+06 );
  $q4 = array(
    2.690530175870899333379843E+03, 
    6.393885654300092398984238E+05, 
    4.135599930241388052042842E+07, 
    1.120872109616147941376570E+09, 
    1.488613728678813811542398E+010, 
    1.016803586272438228077304E+011, 
    3.417476345507377132798597E+011, 
    4.463158187419713286462081E+011 );
  $sqrtpi = 0.9189385332046727417803297;
  $xbig = 4.08E+36;
  
  if ($x<=0 || $xbig < $x) {
	  return 1e10;
  }
  if ($x <= 1e-10) {
	  $res = -log($x);
  } else if ($x<=1.5) {
	  if ($x<$pnt68) {
		  $corr = -log($x);
		  $xm1 = $x;
	  } else {
		  $corr = 0.0;
		  $xm1 = ($x-0.5) - 0.5;
	  }
	  if ($x<0.5 || $pnt68 <= $x) {
		  $xden = 1.0;
		  $xnum = 0.0;
		  for ($i=0; $i<8; $i++) {
			  $xnum = $xnum*$xm1 + $p1[$i];
			  $xden = $xden*$xm1 + $q1[$i];
		  }
		  $res = $corr + ($xm1*($d1 + $xm1*($xnum/$xden)));
	  } else {
		  $xm2 = ($x-0.5)-0.5;
		  $xden = 1.0;
		  $xnum = 0.0;
		  for ($i=0; $i<8; $i++) {
			  $xnum = $xnum*$xm2 + $p2[$i];
			  $xden = $xden*$xm2 + $q2[$i];
		  }
		  $res = $corr + ($xm2*($d2 + $xm2*($xnum/$xden)));
	  }
  } else if ($x<=4.0) {
	  $xm2 = $x-2.0;
	  $xden = 1.0;
	  $xnum = 0.0;
	  for ($i=0; $i<8; $i++) {
		  $xnum = $xnum*$xm2 + $p2[$i];
		  $xden = $xden*$xm2 + $q2[$i];
	  }
	  $res = $xm2*($d2 + $xm2*($xnum/$xden));
  } else if ($x <= 12.0) {
	  $xm4 = $x - 4.0;
	  $xden = -1.0;
	  $xnum = 0.0;
	  for ($i=0; $i<8; $i++) {
		  $xnum = $xnum*$xm4 + $p4[$i];
		  $xden = $xden*$xm4 + $q4[$i];
	  }
	  $res = $d4 + $xm4*($xnum/$xden);
  } else {
	  $res = 0.0;
	  if ($x<= $frtbig) {
		  $res = $c[6];
		  $xsq = $x*$x;
		  for ($i=0;$i<6;$i++) {
			  $res = $res/$xsq + $c[$i];
		  }
	  }
	  $res = $res/$x;
	  $corr = log($x);
	  $res = $res + $sqrtpi - 0.5*$corr;
	  $res = $res + $x*($corr - 1.0);
  }
  return $res;
}

//fcdf(f,df1,df2)
//Returns the area to right of the F-value f for the f-distribution
//with df1 and df2 degrees of freedom (techinically it's 1-CDF)
//Algorithm is accurate to approximately 4-5 decimals
function fcdf($x,$df1,$df2) {
	$p1 = fcall(fspin($x,$df1,$df2));
	return $p1;
}

function fcall($x) {
	if ($x>=0) {
		$x += 0.0000005;
	} else {
		$x -= 0.0000005;
	}
	return $x;
}
function fspin($f,$df1,$df2) {
	$pj2 = M_PI/2;
	$pj4 = M_PI/4;
	$px2 = 2*M_PI;
	$exx = 1.10517091807564;
	$dgr = 180/M_PI;
	
	$x = $df2/($df1*$f+$df2);
	if (($df1%2)==0) {
		return (LJspin(1-$x,$df2,$df1+$df2-4,$df2-2)*pow($x,$df2/2));
	} else if (($df2%2)==0) {
		return (1-LJspin($x,$df1,$df1+$df2-4,$df1-2)*pow(1-$x,$df1/2));
	}
	$tan = atan(sqrt($df1*$f/$df2));
	$a = $tan/$pj2;
	$sat = sin($tan);
	$cot = cos($tan);
	if ($df2>1) {
		$a = $a+$sat*$cot*LJspin($cot*$cot,2,$df2-3,-1)/$pj2;
	}
	if ($df1==1) { 
		return (1-$a);
	}
	$c = 4*LJspin($sat*$sat,$df2+1,$df1+$df2-4,$df2-2)*$sat*pow($cot,$df2)/M_PI;
	if ($df2==1) {
		return (1-$a+$c/2);
	}
	$k = 2;
	while ($k<=($df2-1)/2) {
		$c = $c*$k/($k-0.5);
		$k++;
	}
	return (1-$a+$c);
}

function LJspin($q,$i,$j,$b) {
	$zz = 1;
	$z = $zz;
	$k = $i;
	while ($k<=$j) {
		$zz = $zz*$q*$k/($k-$b);
		$z = $z+$zz;
		$k += 2;
	}
	return $z;
}

//invfcdf(p,df1,df2)
//Computes the f-value with probability of p to the right
//with degrees of freedom df1 and df2
//Algorithm is accurate to approximately 2-4 decimal places
//Less accurate for smaller p-values
function invfcdf($p,$df1,$df2) {
	$v = 0.5;
	$dv = 0.5;
	$f = 0;
	$cnt = 0;
	while ($dv>1e-10) {
		$f = 1/$v-1;
		$dv = $dv/2;
		if (fcdf($f,$df1,$df2)>$p) {
			$v = $v-$dv;
		} else {
			$v = $v+$dv;
		}
		$cnt++;
	}
	echo $cnt;
	return $f;
	
}
?>

